Calibración de probabilidades

Estimación de default mediante modelos de machine learning

Karina Bartolomé

Lic. En economía (UNLP), Esp. Métodos Cuantitativos (UBA)

Organizadora: Natalia Salaberry
CIMBAGE (IADCOM) - Facultad Ciencias Económicas (UBA)

2024-05-13

0.1 Consideraciones previas

  • Se realizará una breve introducción a modelos de aprendizaje automático, pero se recomienda tener algún conocimiento previo para seguir el taller.

  • Durante el seminario se utilizarán los siguientes paquetes (python):

Imports
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, to_rgb
from great_tables import GT, from_column, style, loc
from collections import Counter

# from sklearn import datasets
# Modelado
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.compose import make_column_selector
from sklearn.preprocessing import StandardScaler, MinMaxScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier, plot_tree
from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier
from sklearn.svm import SVC
from sklearn.calibration import calibration_curve, CalibratedClassifierCV
from sklearn.metrics import (
    accuracy_score, precision_score, recall_score, f1_score,
    roc_auc_score,
    roc_curve,
    RocCurveDisplay,
    log_loss,
    recall_score,
    brier_score_loss,
    confusion_matrix,
    ConfusionMatrixDisplay,
)
from IPython.display import display, Markdown

# Funciones adicionales: 
from custom_functions import (
    plot_distribution,
    plot_calibration_models,
    plot_calibration,
    calibrate_model,
    calibration_table,
    calculate_clf_metrics,
)

import sklearn
sklearn.set_config(transform_output="pandas")

Procesamiento de datos:

📦 pandas==2.2.1

📦 numpy==1.26.4

Modelado:

📦 scikit-learn==1.3.0

Visualización y tablas:

📦 matplotlib==3.8.0

📦 seaborn==0.12.2

📦 great_tables==0.5.0

1 Introducción a machine learning

1.1 Machine Learning

En general al hablar de machine learning se hace referencia a diversos tipos de modelos. En la Figura 1 se muestran los distintos tipos de modelos de aprendizaje automático.

En esta presentación se analizará el caso particular de los modelos de clasificación, haciendo énfasis en la Estimación de probabilidad.

%%{init: {'flowchart' : {'curve' : 'linear'}}}%%

flowchart LR
    ml[Machine Learning]:::greenclass1
    supervised[Aprendizaje <br> supervisado]:::greenclass1
    unsupervised[Aprendizaje <br> no supervisado]
    reinforce[Aprendizaje <br> por refuerzo]

    ml_clf[Modelos de <br>clasificación]:::greenclass1
    ml_reg[Modelos de <br>regresión]
    ml_clfclass[Predicción de<br>clase]
    ml_clfrank[Predicción de<br>ordenamiento]
    ml_clfprob[Predicción de <br>probabilidad]:::greenclass2

    ml---temp1[ ]

    temp1-->supervised
    temp1-->unsupervised
    temp1-->reinforce

    supervised---temp2[ ]

    temp2-->ml_clf
    temp2-->ml_reg

    ml_clf---temp3[ ]

    temp3-->ml_clfclass
    temp3-->ml_clfrank
    temp3-->ml_clfprob

    classDef greenclass2 fill:#255255,stroke:#333,stroke-width:2px,color:white;
    classDef greenclass1 fill:#BDCBCC,stroke:#333,stroke-width:2px,color:black;
    style temp1 fill:#FFFFFF00,stroke:#FFFFFF00,height:0px,width:0px;
    style temp2 fill:#FFFFFF00,stroke:#FFFFFF00,height:0px,width:0px;
    style temp3 fill:#FFFFFF00,stroke:#FFFFFF00,height:0px,width:0px;
Figura 1: Tipos de modelos de aprendizaje automático

1.2 ¿Qué es un modelo de clasificación binaria?

Se busca predecir la ocurrencia de un evento a partir de ciertas características observables:

\(P(y=1) = f(X)\)

y: variable que puede tomar 2 valores: 0 o 1

X: matriz de nxm, siendo n la cantidad de observaciones y m la cantidad de variables (o atributos)


Ejemplos populares de modelos de clasificación:

▪️Iris: Clasificación de especies de plantas

▪️Titanic: Clasificación de individuos en sobrevivientes y no sobrevivientes

▪️German Credit: Clasificación de individuos en riesgosos o no riesgosos

2 Caso: German Credit

2.1 Datos

Código
PATH_DATA = 'https://raw.githubusercontent.com/ayseceyda/german-credit-gini-analysis/master/gc.csv'
TARGET = "Risk"
# pd.read_csv(PATH_DATA).to_csv('df_german_credit.csv', index=False)
df = (pd.read_csv(PATH_DATA)
  .drop('Unnamed: 0', axis=1)
  .assign(Risk = lambda x: np.where(x['Risk']=='good',0,1))
)

Se cuenta con un dataset de 1000 observaciones y 10 variables.

Código
def map_color(data):
    return (data[TARGET] == 1).map(
        {True: color_verde_claro, False: 'white'}
    )

(GT(df.sample(6, random_state=123)
        .set_index(TARGET)
        .reset_index()
    )
    .fmt_currency(columns="Credit amount")
    .tab_style(
        style=style.fill(color=map_color), locations=loc.body(columns=df.columns.tolist()),
    )
    .tab_style(
        style=style.text(color='red', weight = "bold"), locations=loc.body(TARGET),
    )
    .tab_options(
        column_labels_background_color=color_verde,
        table_font_names="Times New Roman"
    )
)
Tabla 1: Datos de German Credit (muestra de 6 observaciones)
Risk Age Sex Job Housing Saving accounts Checking account Credit amount Duration Purpose
1 29 male 2 own little little $6,887.00 36 education
1 21 male 2 rent little little $902.00 12 education
0 29 male 1 own moderate $2,333.00 24 furniture/equipment
1 20 female 2 rent little little $2,039.00 18 furniture/equipment
0 35 male 2 own moderate $2,728.00 15 radio/TV
0 23 male 2 own moderate $1,444.00 15 radio/TV

La variable objetivo (target) es Risk, donde el porcentaje de observaciones de clase 1 (riesgosos) es 30.0%

2.2 Particiones

Partición en dataset de entrenamiento, validación y evaluación.

  • Dataset de entrenamiento –> Ajuste del modelo
  • Dataset de validación –> Calibración del modelo
  • Dataset de evaluación –> Métricas del modelo calibrado
Código
y = df[TARGET]
X = df.drop([TARGET], axis=1)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.5, shuffle=True, stratify=y, random_state=42
)

# Lo ideal es utilizar la partición de validación para calibración
X_valid, X_test, y_valid, y_test = train_test_split(
    X_test, y_test, test_size=0.5, shuffle=True, stratify=y_test, random_state=42
)

print(f"N observaciones en entrenamiento: {X_train.shape[0]}")
print(f"N observaciones en validación: {X_valid.shape[0]}")
print(f"N observaciones en evaluación: {X_test.shape[0]}")
N observaciones en entrenamiento: 500
N observaciones en validación: 250
N observaciones en evaluación: 250

2.3 Modelo de clasificación simple

Se consideran 2 variables para ajustar un modelo de arbol de decisión con máxima profundidad=2

\[ P(\text{Risk}=1) = f(\text{Credit amount}, \text{Age}) \]

Código
subset_cols = ['Age', 'Credit amount']

(GT(pd.concat([y_train, X_train], axis=1)
        .sample(6, random_state=123)
        .set_index(TARGET)
        [subset_cols]
        .reset_index()
    )
    .fmt_currency(columns="Credit amount")
    .tab_style(
        style=style.fill(color=map_color), locations=loc.body(columns=df.columns.tolist()),
    )
    .tab_style(
        style=style.text(color='red', weight = "bold"), locations=loc.body(TARGET),
    )
    .tab_options(
        column_labels_background_color=color_verde,
        table_font_names="Times New Roman"
    )
)
Tabla 2: Datos de German Credit (2 variables)
Risk Age Credit amount
0 37 $1,287.00
1 27 $8,648.00
1 34 $1,337.00
0 48 $2,134.00
1 48 $1,082.00
0 53 $2,424.00

Ajuste del modelo:

Código
X_train_subset = X_train[subset_cols].copy()
clf = DecisionTreeClassifier(max_depth=2).fit(X_train_subset, y_train)
clf
DecisionTreeClassifier(max_depth=2)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
Código
colors = [color_verde, 'red']
labels = ['No riesgoso','Riesgoso']
plt.figure(figsize=(6,6))
arbol = plot_tree(
    clf,
    filled=True,
    impurity=False,
    feature_names=X_train_subset.columns.tolist(),
    class_names=labels,
    fontsize=8
)
for i, impurity, value in zip(arbol, clf.tree_.impurity, clf.tree_.value):
    # let the max value decide the color; whiten the color depending on impurity (gini)
    r, g, b = to_rgb(colors[np.argmax(value)])
    f = impurity * 2 # for N colors: f = impurity * N/(N-1) if N>1 else 0
    i.get_bbox_patch().set_facecolor((f + (1-f)*r, f + (1-f)*g, f + (1-f)*b))
    i.get_bbox_patch().set_edgecolor('black')
plt.show()
Figura 2: Modelo de árbol de decisión

2.4 Evaluación de un modelo de clasificación

Se utiliza el método clf.predict()para predecir la clase (0,1) según el modelo (clf). A partir de la clase predicha se obtiene la matriz de confusión:

Ver código
preds = pd.DataFrame({'Valor observado':y_test,'Valor predicho':clf.predict(X_test[subset_cols])})
display(GT(preds
    .groupby(['Valor observado','Valor predicho'], as_index=False).size()
    .pivot(index='Valor observado', columns='Valor predicho', values='size')
    .rename({0:'0',1:'1'}, axis=1)
    .reset_index()
    )
    .tab_spanner('Valor predicho', columns = ['0','1'])
    .data_color(
        columns=['0','1'],
        domain=[0,X.shape[0]],
        palette=[color_verde_claro, 'red'],
        na_color="white",
    )
)
Tabla 3: Matriz de confusión
Valor observado Valor predicho
0 1
0 172 3
1 71 4

Además de predecir una clase, sklearn permite predecir una “probabilidad” con clf.predict_proba():

Ver código
y_pred_proba = clf.predict_proba(X_test[subset_cols])[:, 1]
preds = pd.DataFrame({
    'y_obs':y_test,
    'predict_proba':y_pred_proba
})
GT(preds.sample(2, random_state=42)).fmt_number('predict_proba')
Tabla 4: clf.predict_proba(), muestra de 2 observaciones
y_obs predict_proba
0 0.25
1 0.33

ROC AUC = 0.64, Log loss = 0.59, Brier loss = 0.2
Ver código
distrib=plot_distribution(data=preds, pred_column='predict_proba', class_0_color=color_verde)
distrib
Figura 3: Distribución de “probabilidad” predicha

2.5 Evaluación de un modelo de clasificación (Cont.)

En la Figura 4 se presentan distintas métricas de clasificación. En este caso el foco estará sobre las métricas basadas en probabilidad.

%%{init: {'flowchart' : {'curve' : 'linear'}}}%%

flowchart LR
    ml[Modelo de <br>clasificación]
    metrics[Métricas]

    ml_clfclass[Predicción de<br>clase]
    ml_clfrank[Predicción de<br>ordenamiento]
    ml_clfprob[Predicción de <br>probabilidad]:::greenclass2
    
    ml-->metrics
    metrics---temp1[ ]

    temp1-->ml_clfclass
    temp1-->ml_clfrank
    temp1-->ml_clfprob

    ml_clfclass---temp2[ ]
    temp2-->accuracy
    temp2-->precision
    temp2-->etc

    ml_clfrank---temp3[ ]
    temp3-->roc_auc
    
    ml_clfprob---temp4[ ]
    temp4-->log_loss:::greenclass2
    temp4-->brier_loss:::greenclass2

    classDef greenclass2 fill:#255255,stroke:#333,stroke-width:2px,color:white;
    classDef greenclass1 fill:#BDCBCC,stroke:#333,stroke-width:2px,color:black;
    style temp1 fill:#FFFFFF00,stroke:#FFFFFF00,height:0px,width:0px;
    style temp2 fill:#FFFFFF00,stroke:#FFFFFF00,height:0px,width:0px;
    style temp3 fill:#FFFFFF00,stroke:#FFFFFF00,height:0px,width:0px;
    style temp4 fill:#FFFFFF00,stroke:#FFFFFF00,height:0px,width:0px;
Figura 4: Métricas de modelos de clasificación

2.6 Evaluación de un modelo de clasificación (Cont.)

Siendo \(N\) la cantidad de observaciones, \(y_i\) el valor observado para la obseravción \(i\) y \(\hat{y_i}\) el valor predicho para la observación \(i\), se definen las funciones de pérdida:

3 Comparación de modelos de ML

3.1 Preprocesamiento

Muchos modelos de aprendizaje automático requieren que los datos sean numéricos. Para ello, se construye un pipeline de preprocesamiento de variables, diferenciando el procesamiento según tipo de variables: categóricas y numéricas.

🔗 sklearn.preprocessing

Ver código
numeric_transformer = Pipeline([
    ('impute', SimpleImputer(strategy='median')),
    ('scaler', MinMaxScaler())
])

categorical_transformer = Pipeline([
    ('ohe',OneHotEncoder(
        min_frequency=0.05,
        handle_unknown='infrequent_if_exist',
        sparse_output=False)
    )
])

preproc = ColumnTransformer([
    ('num', numeric_transformer,
      make_column_selector(dtype_include=['float','int'])),
    ('cat', categorical_transformer,
      make_column_selector(dtype_include=['object','category']))
], verbose_feature_names_out=False)
ColumnTransformer(transformers=[('num',
                                 Pipeline(steps=[('impute',
                                                  SimpleImputer(strategy='median')),
                                                 ('scaler', MinMaxScaler())]),
                                 <sklearn.compose._column_transformer.make_column_selector object at 0x146486260>),
                                ('cat',
                                 Pipeline(steps=[('ohe',
                                                  OneHotEncoder(handle_unknown='infrequent_if_exist',
                                                                min_frequency=0.05,
                                                                sparse_output=False))]),
                                 <sklearn.compose._column_transformer.make_column_selector object at 0x146484580>)],
                  verbose_feature_names_out=False)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

En la Tabla 5 se visualizan los datos luego del preprocesamiento. Esta es la matriz que será utilizada para ajustar múltiples modelos de clasificación.

Código
display(GT(preproc.fit_transform(X_train).sample(2).round(2))
    .tab_options(
        column_labels_background_color=color_verde,
        table_font_names="Times New Roman"
    )
)
Tabla 5: Datos de German Credit preprocesados (muestra de 2 observaciones)
Age Job Credit amount Duration Sex_female Sex_male Housing_free Housing_own Housing_rent Saving accounts_little Saving accounts_moderate Saving accounts_quite rich Saving accounts_rich Saving accounts_nan Checking account_little Checking account_moderate Checking account_rich Checking account_nan Purpose_business Purpose_car Purpose_furniture/equipment Purpose_radio/TV Purpose_infrequent_sklearn
0.07 0.67 0.49 0.57 1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0
0.48 0.67 0.43 0.79 0.0 1.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0

3.2 Modelado

Se ajustan distintos tipos de modelos reutilizando el mismo pipeline de preprocesamiento de datos:

3.3 Evaluación

Se calculan las métricas de cada uno de los modelos:

Tabla 6: Métricas de los modelos en la partición de evaluación
Modelo Accuracy Precision Recall F1 ROC AUC Log Loss Brier Loss
Hist gradient boosting 0.72 0.53 0.53 0.53 0.74 0.93 0.22
Decision tree 0.69 0.49 0.47 0.48 0.63 11.1 0.31
Logistc regression 0.74 0.59 0.45 0.51 0.75 0.53 0.18
SVC 0.73 0.6 0.33 0.43 0.72 0.55 0.19
Logistc regression (balanced) 0.64 0.44 0.69 0.54 0.75 0.6 0.21
Hist gradient boosting (balanced) 0.71 0.52 0.59 0.55 0.74 0.98 0.23

3.4 Problema organizacional

Se utiliza el modelo para generar predicciones sobre las observaciones de test (Tabla 7), agrupando las predicciones en 5 intervalos de probabilidad. La Tabla 8 contiene la cantidad de observaciones por bin y el promedio de y_obs e y_pred. En la Figura 6 se visualiza la probabilidad promedio por bin en el eje x y la fracción de clase positiva observada en el eje y.

Ver código
preds_df = (pd.DataFrame({
        'y_obs':y_test,
        'y_pred':pipe_hist.predict_proba(X_test)[:,1]
    })
    .assign(bin=lambda x: pd.qcut(x['y_pred'],q=5, precision=5))
)
preds_df['bin'] = pd.IntervalIndex(preds_df['bin']).map(lambda x: pd.Interval(round(x.left, 2), round(x.right, 2)))
(GT(preds_df
    .sample(3, random_state=123).round(4))
    .tab_options(
        column_labels_background_color=color_verde,
        table_font_names="Times New Roman"
    )
)
Tabla 7: Predicciones
y_obs y_pred bin
0 0.1315 (0.02, 0.23]
0 0.0 (-0.0, 0.0]
1 0.6428 (0.23, 0.84]
Ver código
(GT(preds_df.groupby('bin')
    .agg(
        N = ('y_obs','count'), 
        frac_positive=('y_obs','mean'),
        avg_pred = ('y_pred','mean')
    ).round(2).reset_index()
    )
    .tab_options(
        column_labels_background_color=color_verde,
        table_font_names="Times New Roman"
    )
)
Tabla 8: Bins de las predicciones
bin N frac_positive avg_pred
(-0.0, 0.0] 50 0.08 0.0
(0.0, 0.02] 50 0.14 0.01
(0.02, 0.23] 50 0.28 0.08
(0.23, 0.84] 50 0.4 0.51
(0.84, 1.0] 50 0.6 0.95
Ver código
from sklearn.calibration import calibration_curve
prob_true, prob_pred = calibration_curve(
    y_test, pipe_hist.predict_proba(X_test)[:,1], 
    n_bins=5,
    strategy='uniform'
)
plt.figure(figsize=(3,3))
sns.lineplot(x=prob_pred, y=prob_true, color='red', label='Modelo')
sns.scatterplot(x=prob_pred, y=prob_true, color='red')
plt.axline([0, 0], [1, 1], c=color_verde)
plt.legend()
plt.xlabel("Probabilidad predicha")
plt.ylabel("Fracción de positivos");
Figura 6: Gráfico de calibración (Reliability diagram)

4 Calibración de probabilidades

4.1 Calibración de probabilidades

Objetivo:

Se busca encontrar una función que ajuste la relación entre los scores predichos por el modelo (predict_proba) y las probabilidades reales

\(P(y_i=1)= f(z)\)

Siendo: \(P(y_{i}=1)\) la probabilidad calibrada para el individuo \(i\) y \(z_{i}\) el output del modelo no calibrado (score)

Un clasificador binario bien calibrado debería clasificar de forma tal que para las observaciones que predice un score (predict_proba) = 0.8, se espera que un 80% de las observaciones correspondan a la clase positiva (1).

Referencias:

4.2 Calibración de probabilidades (Cont.)

Métodos de calibración

Calibración Sigmoide (Platt scaling)

Calibración Isotónica

Calibración Beta (Nuevo)

Calibración Spline (Nuevo)

4.3 Calibración sigmoide (Platt scaling)

Este método asume una relación logística entre los scores (z) y la probabilidad real (p):

\(log(\frac{p}{1-p})=\alpha+\beta(z_{i})\)

\(P(y_{i}=1) = \frac{1}{1+(e^{-(α+β(z_{i}))})}\)

Siendo: \(y_{i}\) el valor observado para el individuo \(i\) y \(z_{i}\) el output del modelo no calibrado


Se estiman 2 parámetros (\(\alpha\) y \(\beta\)), como en una regresión logística.

  • Requiere pocos datos

  • Es útil cuando el modelo no calibrado tiene errores similares en predicción de valores bajos y altos


Referencias:

  • Platt, John. (2000). Probabilistic Outputs for Support Vector Machines and Comparisons to Regularized Likelihood Methods. Adv. Large Margin Classif, Volume 10 (pp. 61-74).

4.4 Calibración sigmoide (Platt scaling) (Cont.)

Ver código
pipe_sigmoid_hist, _, preds_hist = calibrate_model(
    pipe=pipe_hist,
    X_valid=X_valid, y_valid=y_valid, 
    X_test=X_test, y_test=y_test,
    cal_sigmoid=True,
    cal_isotonic=False,
    model_name=''
)
plot_calibration(
    preds=preds_hist,
    y_pred_prob_calibrated='pred_sigmoid',
    model_name=''
)

4.5 Calibración isotónica

Se ajusta un regresor isotónico no paramétrico, que produce una función escalonada no decreciente:

\(\sum_{i=1}^{n}(y_i - \hat{f_i})^2\)

Sujeto a \(\hat{f_i}\) >= \(\hat{f_j}\) siempre \(f_i\) > \(f_j\) \(y_i\): etiqueta verdadera de la obseravción \(i\) y sea la salida del clasificador calibrado para la muestra (es decir, la probabilidad calibrada)

4.6 Calibración isotónica (Cont.)

4.7 Tabla de calibración

4.8 Comparativa de métricas

5 Comentarios finales

5.1 Comentarios finales

Se busca un modelo en donde la probabilidad promedio de cada bin se corresponda con el % de clase positiva observado en ese bin según las predicciones del modelo.

Esto es útil para la toma de decisiones, ya que permite establecer punto de cortes diferenciales en función de la probabilidad.

5.2 Referencias / Recursos

Probability Calibration Workshop - PyData 2020

5.3 Contacto

karinabartolome

@karbartolome

@karbartolome

Blog